classdef QA  <QuesNewtonBasicMethod
    properties 
        fmax
        Indx
    end
    methods       
          
        function [Solution obj] = Run(obj)
            %% Culot Algorithm ;
            %% A quasi-Newton algorithm for first-orderSaddle Points
            %% PDF: -> 'D:\Saddle Points Literature\Transition States\Transition
            %% State\Deterministic Methods\Single ended\A quasi-Newton algorithm for first-orderSaddle Points.pdf'

            obj.Dims = 2;
            Solution = [];
            %X = [0; 0];
            global X;
            %X = [-1.5; -0.5];
            %% 1) Calculate the energy and the derivatives at the initial geometry
            [obj.E, obj.Grad, obj.H]= obj.CalculateRealEnergyDerivatives(X);
            isSolution = false;
            obj.Rk = 0.2;
            
            %% 2) Diagnolize the hessian matrix and choose an eigenvector
            %% vj
             if isnan(det(obj.H)) || isinf(det(obj.H)) || det(obj.H) ==0
                    fprintf('\n Singular Matrix encountered during the search\n');
                    return;
             end
            
            [obj.eigVectors eigValsM] = eig(obj.H);
            
            obj.eigVals = zeros(1,obj.Dims)
            for i = 1: obj.Dims
               obj.eigVals(i) =  eigValsM(i,i);
            end    
            
            DI = DirectionIdentification;
            [minEigVal obj.Indx]= DI.FindMinIndex(obj.eigVals, obj.Dims);
            maxEigVal = max(obj.eigVals);
            obj.counter = 0;
            for i=1: obj.MaxIterations                             
                %% 3) Calculate Step Vector 
                [obj.n_pos, obj.n_zero, obj.n_neg] = obj.DetermineInertiaOfHessian(...
                    obj.eigVals, obj.Dims);
                obj.counter = obj.counter + 1
                F = [0, maxEigVal, -minEigVal];
                obj.fmax = max(F);
                
                Step = obj.MoveAStep();
                
                %% 5) Calculate the quadratic energy function approximation    
                %% 6) Calculate the energy at the point X + D upadte the trust radius
                E1 = obj.CalculateQuadraticEnergy(Step, obj.E, obj.Grad, obj.H);    
                
                Ce = obj.CalculatePredictedChangeInEnergy( obj.Grad, Step, obj.H);
                TrstRedius = TrustRegion;
                Rk = TrstRedius.UpdateTrustRaduis(obj.Rk, E1, obj.E, Ce, Step, obj.Grad, obj.H);                             
              
              
                %% is the current configruation is transition states or first order
                %% saddle points
                %% 7.2.1) Yes: report result
                
                GradRMS = norm(obj.Grad);
                MaxGrad = max(obj.Grad);
                MaxStep = max(Step);
                StepRMS = norm(Step);
                
                if StepRMS == Inf
                    continue;
                end

                isSolution = obj.isATransitionState(GradRMS, MaxGrad,...
                    MaxStep, StepRMS);
                if isSolution                     
                    break;
                end

                %% 7.2.2) No:
                %% Update the geometry X, Calculate the gradient vector at the new
                %% geometry, update the hessian matrix     
                X = X + Step;
                obj.E = E1;   
                Gk = obj.Grad;
                obj.Grad = Gradient(X);
                Gk1 = obj.Grad;
                obj.H = obj.UpdateHessian(Step, Gk1, Gk,obj.H);                
                OldVector = obj.eigVectors(:,obj.Indx);                
                if isnan(det(obj.H)) || isinf(det(obj.H)) || det(obj.H) ==0
                    fprintf('\n Singular Matrix encountered during the search\n');
                    return;
                end
                [obj.eigVectors eigValsM] = eig(obj.H);         
                obj.eigVals = zeros(1,obj.Dims);
                for i = 1: obj.Dims
                   obj.eigVals(i) =  eigValsM(i,i);
                end
                minEigVal = DI.FindMinIndex(obj.eigVals, obj.Dims);
                maxEigVal = max(obj.eigVals);
                obj.Indx = DI.IdentifyDirection( OldVector, obj.eigVectors, obj.Dims);
            end
            if isSolution
                Solution = X
            else
                Solution = [];
            end
            
        end
        
        function Step = MoveAStep(obj) 
            Step = obj.CalculateQausiStep();
%             if norm(Step) > obj.Rk
%                 Step  = Step * obj.Scaler;                
%             end            
            if obj.n_neg ~= 1 || norm(Step) > obj.Rk
                fprintf('\n Augemented QuasiNetwon Step ');                
                Step = obj.CalculateStep();
            end
            
        end
        
        
         function Step = CalculateStep(obj)            
            y = obj.IdentifyLamba();
            Step = zeros(size(obj.Grad));
            for i = 1: obj.Dims                
                Grad_ = obj.eigVectors(:,i) * obj.Grad';
                if i == obj.Indx
                    Grad_ = Grad_/(obj.eigVals(i)- y ) ;
                else 
                    Grad_ = Grad_/(obj.eigVals(i)+ y);
                end
                Step = Step - Grad_ * obj.eigVectors(:,i);    
            end
         end
       
        function lamba = IdentifyLamba(obj)
            NR = NewtonRaphsonMethod ;
            Func =@ obj.EvaluateQALamba;
            d = sort(obj.eigVals);
            InitialGuess = 0.25*(d(1) + d(2)/2);                         
            %% we need to check it without constraints
            lamba = NR.Run(obj.NWMaxIterations, ...
                 InitialGuess, Func, obj.Dims);   
        end        
        
        
       function Fit = EvaluateQALamba(obj,y)
             Step = zeros(size(obj.Grad));
             for i = 1: obj.Dims                
                Grad_ = obj.eigVectors(:,i) * obj.Grad';
                if i == obj.Indx
                    Grad_ = Grad_/(obj.eigVals(i)- y ) ;
                else
                    Grad_ = Grad_/(obj.eigVals(i)+ y);
                end
                Step = Step - Grad_' * obj.eigVectors(:,i);
            end
            
            FinalFunc = 0;
            for i = 1: obj.Dims                
                FinalFunc = FinalFunc + Step(i)^2;
            end
            FinalFunc = sqrt(FinalFunc);
            Fit = (FinalFunc - obj.Rk)^2;
        end   

    end
end
    


